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The past decade has seen a revived interest in the unavoidable or intrinsic noise in biochemical 
and genetic networks arising from the finite copy number of the participating species. That is, 
rather than modeling regulatory networks in terms of the deterministic dynamics of concentrations, 
we model the dynamics of the probability of a given copy number of the reactants in single cells. 
Most of the modeling activity of the last decade has centered on stochastic simulation of individual 
realizations, i.e., Monte-Carlo methods for generating stochastic time series. Here we review the 

r- h mathematical description in terms of probability distributions, introducing the relevant derivations 

and illustrating several cases for which analytic progress can be made either instead of or before 

^Nj turning to numerical computation. 
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I. INTRODUCTION 

Cellular processes rely on biochemical signaling, i.e. chemical interactions among individual molecules. Theoretical 
descriptions of biochemical processes rely on considering the concentration changes of the molecules involved in the 
reactions, and specifying the types of regulatory interactions between them. The theoretical tools used to describe 
many types of biochemical reactions share many similarities. In this chapter we give an overview of the analytical 
approaches to study biochemical kinetics using the example of small gene regulatory networks. 

The regulation of genes by transcription factor proteins is an intrinsically stochastic process, owing to the small 
numbers of copies of molecules involved. With the development of imaging techniques in molecular biology, we are 
able to observe directly the fluctuations in the concentrations of proteins and mRNAs and, by measuring the intensity 
profiles of fluorescence markers, measure full probability distributions [1-4]. Experiments over the last decade have 
shown that in fact gene regulation is a noisy process, and noise can propagate in gene networks. 

Many methods for solving the resulting stochastic equations rely on computer simulations. The efficiency of these 
methods has been greatly advanced in the last several years [5-16] . However numerical simulations are naturally 
limited to a specific choice of parameters, and changing the parameters requires a completely new calculation. Fur- 
thermore they suffer from the curse of dimensionality: the computational runtime grows prohibitively as the number 
of species increases These problems can be bypassed by developing analytical approaches, which often require certain 
approximations. 

This chapter is intended as a tutorial on theoretical descriptions of biochemical kinetics. For clarity of exposition, 
we start by introducing the simplified kinetic description of the production of proteins. At the level of molecular 
kinetics, this is the simplest type of process that can occur. We first consider the deterministic description of the 
system and then introduce noise. After familiarizing the reader with different levels of description, we discuss models 
of regulation. We also present a wide spectrum of analytic tools used to find the steady state probability distributions 
of protein copy numbers in the cell. We point out that while for concreteness we focus on the case of gene regulation, 
the methods presented in this chapter are very general. 

The results presented in this review are not new; some can be found in textbooks [17-19], while others have been 
derived more recently in the context of gene regulation [4, 20-38]. The goal of this review is to give the reader who is 
not familiar with analytical methods for modeling stochastic gene regulation an overview of the mathematical tools 
used in the field. Naturally, we are not able to cover all the developments in the field, but we hope to give the 
reader a useful starting point. For concreteness we also limit our discussion to the case of small gene networks and 
do not discuss approximations used to describe larger networks, which is currently an active area of research in many 
communities [10, 11, 15, 16, 28, 29]. 



II. SIMPLE BIRTH-DEATH PROCESS 

In this section we describe a simple birth death process for one species, absent of regulation. We first remind the 
reader of the deterministic description given by chemical kinetics. Next we introduce the full probabilistic description 
(given by the master equation), and, restricting our discussion to the simplest case, show how the deterministic 
equations arise as the dynamics of the mean. After calculating the variance and, we introduce the generating function 
formalism, useful in solving the master equation. This formalism may also be described using raising and lowering 
operators familiar to physicists from quantum mechanics. Finally we relate this description to the Fokkcr-Plank 
equation, the analogous master equation for continuous state variables (e.g., real-valued coordinates rather than the 
integer- valued copy number). These results will be used in later sections when we introduce autoregulation and 
regulation among different species. 

A. Deterministic description: the kinetic rate equation 

In the simplest case, the number of copies n of a protein species X can change either due to the production of a 
protein, which occurs at a constant rate g, or due to degradation of a protein, which occurs at a constant rate r: 

0^X (1) 

r 

Here we condense the complicated molecular machinery of transcription, translation, and protein modification into a 
single constant production rate. Similarly we do not specify the molecular processes which have led to degradation; 
for simplicity we imagine either dilution due to cell division or active degradation with a constant rate. 

When n is large, the dynamics for the mean (n) are well approximated by the continuous, deterministic description 
provided by the kinetic rate equation for the concentration c = (n)/V in a volume V: 



The solution of Eqn. 2 is 
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where c(0) is the concentration of proteins at initial time t = 0. In steady state the mean number of proteins is simply 
the ratio of the production and degradation rates, (n) = cV = g/r, which is easily seen by taking either dc/dt = in 
Eqn. 2 or t — > oo in Eqn. 3. 

B. Introducing noise: the master equation 

The general probabilistic description of chemical reactions, respecting the finite copy number of the reactants, is the 
master equation, which specifies the rate of change of p ni the probability that there are n = (ni, ri2, . . . , Ul) copies of 
the L reactants. The macroscopic, deterministic description in terms of chemical kinetics is recovered by considering 
the dynamics for the mean concentration of these reactants, i.e. the vector {n)/V = {{n\)/V, (n2) /V, . . . (til)/V) = 
(ci, C2, . . . Cl,). Such a description in terms of a summary statistic necessarily ignores a great amount of information, 
including all information about fluctuations about these mean values. 

The general form of the master equation (c.f. Appendix A), set by conservation of total probability, is 



Ei 



' (t)Pn> ~ Wnn' (*)Pn] , (4) 



For the case of only one species, we need specify only the dynamics of n = n. Restricting further to the case of a 
simple birth-death process, there are only two nonzero contributions from the transition matrices w nn >, given by 

w n +i,n(i) = g, (5) 

w n -i >n (t) = rn, (6) 



all other values of w nn i being 0. Under this restriction Eqn. 4 reduces to the familiar 

Pn = -(JPn - rnp n + gp n -i + r(n + l)p n +i- (7) 

Qualitatively, the four terms on the right hand side represent how the probability of having n proteins can either 
(i) decrease in time, if there are n initially and one is either produced (first term) or degraded (second term), or 
(ii) increase in time, if there are either n — 1 initially and one is produced (third term) or n + 1 initially and one is 
degraded (fourth term). 

The dynamics of the mean number of proteins (n) = X^^Lo n P n are rea dily obtained (see Appendix B) , as 

din) . . ,„. 

■y-=9-r(n). (8) 

Eqn. 8 shows that the dynamics of the mean of the protein distribution reproduces the kinetic rate equation, Eqn. 2. 
From here on, we will re-scale time t by the degradation rate r, such that rt — > t, and define g = g/r, making the 
master equation 

Pn = -9Pn - np n + gp n -l + (n + l)p n+1 . (9) 

1. Steady state solution 

The master equation can be rewritten in terms of shift operators E + and E~ which increase and decrease the 
number of proteins by one, respectively [17], i.e. 

E + fn = /„+!, (10) 

E-fn = fn-1, (11) 

for any function /„. We begin by writing Eqn. 9 in terms of only E + to make clear that, as is the case for any 
one-dimensional master equation, its steady state can be found iteratively. In terms of E + , 

Pn = (E + - l){np„ - gp n -i). (12) 

Setting p n = for steady state requires that the second term in parentheses vanish, giving for n > the recursive 
relation 

Pn=Pn-i-- (13) 

n 

Starting with n = 1 and computing the first few terms reveals the pattern 

g n 

Pn = — 7P0 (14) 
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where po = e 9 is set by normalization: 

1 = ^2p^=Po^2zJ =Poe 9 . (15) 



.x: 



n=0 n=0 

Thus the steady-state probability of having n proteins is the Poisson distribution, 



Pn = -^\ (16) 

n\ 

with parameter g, the ratio of production to degradation rates. 

We remind the reader that the variance a 1 = (n 2 ) — (n) 2 of the Poisson distribution is equal to its mean (n) 
(Appendix C). Therefore the standard deviation a over the mean falls off like 

— = ^ = — (17) 

(n) (n) ,/M' l ] 



demonstrating that the relative effect of fluctuations diminishes for large protein number. In Appendix D we show 
that in the limit of large protein number the steady state asymptotes to a Gaussian distribution with mean and 
variance equal to g. 

For comparison with the other representations of the master equation described below, we now write Eqn. 9 in 
terms of both shift operators E + and E~ : 

i>n = {E+ -\){n-gE-)p n . (18) 

Eqn. 12 can be rewritten slightly by inserting between the parenthetic terms the unit operator 1 = E~ E + and 
distributing the E~ to the left and the E + to the right, giving 

p n = -(E--l)[(n + l)E + -g]p n . (19) 

where the negative sign has been factored out of the first parenthetic term. 

2. Generating function representation 

Not all master equations are solvable by straightforward iteration. A more generalizable way to solve a master 
equation is by the introduction of a generating function [17]. Here we demonstrate the generating function approach 
on the birth-death master equation, Eqn. 9. 

The generating function G is defined as 



G(x,t) = J2Pn(t)x n , (20) 



n=0 

a power-series expansion in a continuous variable x whose coefficients are the probabilities p n . Since a; is a variable 
we introduce, we note that defining x = e lK makes clear that the generating function is equivalent to the Fourier 
transform of p n in protein number (this point is further developed in Appendix E). The probability distribution may 
be recovered via the inverse transform 

p n (t) = -.d2G(x,t)\ x = . (21) 

nl 

Additionally we note that the ^th moment may be generated by 

(n')=3?GU =1) (22) 

(which is the reason for the generating function's name). 

The utility of the generating function for solving a master equation is that it turns an infinite set of ordinary 
differential equations (e.g. Eqn. 9) into a single partial differential equation. To see this here we multiply Eqn. 9 by 
x n and sum over n to obtain (see Appendix B for the detailed derivation): 

G = -(x-l)(d x -g)G. (23) 

We see immediately from comparison of Eqns. 19 and 23 that the representations of operators E~ and (n + 1)E + in 
x space are x and d x , respectively. 

In steady state (G = 0) Eqn. 23 must satisfy (d x — g)G = 0, which is solved by 

G(x) = G(0)e 9X , (24) 

where G(0) = e~ 9 is set by normalization: 

oo oo 

G(0)e 9 = G(l) = Y, Pn{l) n = E Pn = I- (25) 

n=0 n=0 

The steady state distribution is recovered via inverse transform: 



Pn — . V x 



(■:■' 



= e~° 9 -, (26) 

x=o nl 



a Poisson distribution with parameter g, as before (Eqn. 16). 

The time-dependent solution to Eqn. 23 may be obtained using the method of characteristics, in which one looks 
for characteristic curves along which the partial differential equation becomes an ordinary differential equation. The 
curves are parameterized by s, such that G(x,t) = G(x(s),t(s)) = G(s), and we look for the solution to the ordinary 
differential equation for G(s), having eliminated x and t from the problem. First Eqn. 23 can be rewritten more 
explicitly as 



Using the chain rule on G(s) gives 



dG f)C 
(x-l) 9 G=- + -(x-l). (27) 



(28) 



dG 

ds 


dGdt 
dt ds 


dGdx 
dx ds 


dt 
ds 


= 1, 




dx 
ds 


= x — 1 




dG 
ds 


= (x-1 


)gG. 


ss of g 


enerality) 
y = VaP 1 - 


. Ther 



Consistency of Eqns. 27 and 28 requires 

(29) 

(30) 

(31) 

Eqn. 29 implies s = t (we set to = without loss of generality). Therefore, defining y = x — 1, Eqn. 30 implies 

(32) 
Finally, straightforward integration of Eqn. 31 yields 

G = G exp gy / dt' e*' = G exp [gy (e* - 1)] = F(y ) cxp [gy e*] . (33) 

Jo J 

where we define F(yo) = Goe~ 9yo . We may expand -F(yo) as a function of its argument, i.e. F(yo) — X^=o Aj^o f° r 



some coefficients Aj, such that 



3 

-t _ t„ i \„-t 



Inserting y = ye = (x — l)e (Eqn. 32), we obtain 

G(x,t) = ^ j A j e- jt {x-lYe a{x - 1) . (35) 

3 

For t — > oo only the j = term survives and we recover the steady state generating function (Eqn. 24), 

G(x) = coe 9 ^ 1 ' (36) 

where cq = 1 by normalization. Eqn. 35 constitutes the full time-dependent solution to the birth death process; 
the time-dependent distribution can be retrieved using the inverse transform, Eqn. 21, and the coefficients Aj arc 
computable from the initial distribution p n (0), which will be made explicit in the next section. 

Since Eqn. 23 is linear in G, one may arrive at its solution in a second, elegant way by expanding in the eigenfunctions 
of the x-dependent operator. That is, writing Eqn. 23 as G — ~CG, where 

C = (x-l)(d x -g), (37) 

we expand G in eigenfunctions 4>j(x) with time-dependent expansion coefficients Gj{t), 

G(x > t) = J2 G At)<!>3(x)> (38) 



where the <fij(x) satisfy 

dj>j = Xrfj (39) 

for eigenvalues Xj. Substituting Eqn. 38 into Eqn. 23 gives Y] ■ Gj4>j = — Y] ■ XjGj(f>j, which, by orthogonality of the 
4>j, yields the set of ordinary differential equations Gj = —XjGj, solved by 

G i {t) = A i e~ x > t (40) 

for some constants Aj. Comparing the result, 

G{x,t) = Y t A i e~^ t <l> i {x\ (41) 

i 

with Eqn. 35 reveals the forms of the eigenvalues 

A j= .?g {0,1,2,...} (42) 

and the eigenfunctions 

<f> j {x) = {x-l) j e 9{x - 1 \ (43) 

facts that we will confirm in the next section using operator methods. 

3. Operator representation 

We now introduce an a representation of the master equation in terms of raising and lowering operators. This 
representation makes the solution more elegant, yielding a simple algebra which allows calculation of projections 
between the spaces of protein numbers and eigenfunctions without explicit computation of overlap integrals; it also 
lays the formal groundwork for solving models of multidimensional regulatory networks. The formalism was first used 
for diffusion by Doi [39] and Zel'dovich [40] and later developed by Pcliti [41]. 

As before the generating function is defined as an expansion in a complete basis indexed by protein number n in 
which the expansion coefficients arc the probabilities p n . Within the operator formalism, this basis is represented as 
the set of states In), i.e. 



|G(t)> = X> n (t)|n>. (44) 



n=0 

Here we have adopted state notation commonly used in quantum mechanics [42]; the previous representation (Eqn. 
20) is recovered by projecting onto this equation the conjugate states (x\ with the provisions 

(x\G(t)) = G(x,t), (45) 

(x\n) = x n . (46) 

In Appendix E we show how the orthonormality of the \n) states, 

(n\n') = S rm > ) (47) 

dictates the form of the conjugate state: 

Projecting (n\ onto Eqn. 45 and using Eqn. 47 gives the inverse transform: 

(n\G(t))= Pn (t). (49) 

The equation of motion in the operator representation is obtained by summing the master equation (Eqn. 9) over 
n against \n) (see Appendix B), giving 

\G) = -£\G), (50) 



where 

C = (a + -l)(a--g). (51) 

Just as in the operator treatment of the quantum harmonic oscillator [42] , the operators a + and aT here raise and 
lower the protein number by 1 respectively, i.e. 

a+\n) = |n+l), (52) 

a~\n) = n\n — 1) (53) 

(note, however, that the prefactors here, 1 and n, are different than those conventionally used for the harmonic oscil- 
lator, \Jn + 1 and y/n, respectively). Comparison of Eqns. 19, 23, and 51 makes clear the following correspondences 
among the master equation, generating function, and operator representations respectively: 

E~ o x o a + , (54) 

(n+l)E+ <-> d x «• a~. (55) 

While it might seem strange that the down-shift operator E~ corresponds to the raising operator d + (and the up-shift 
operator to the lowering operator), in Appendix F we show that, just as with the quantum harmonic oscillator, a + 
and a~ lower and raise protein number respectively when operating to the left, i.e. 

(n\a+ = (n-l|, (56) 

(n\a- = (n + l)<n+l|, (57) 

which makes the correspondence more directly apparent. Finally, again as with the quantum harmonic oscillator, the 
raising and lowering operators enjoy the commutation relation (see Appendix F) 

[a-,a+]=i, (58) 

and a + a~ acts as a number operator, i.e. a + ar\n) = n\n). 

Eqn. 51 shows that the full operator C factorizes, suggesting the definition of the shifted operators 

b+ = a + - 1, (59) 

b~ = cT - g. (60) 

Since b + and b~ differ from a + and aT 1 respectively, by scalars, they obey the same commutation relation, i.e. 



b~y 



= 1- (61) 



Since the equation of motion (Eqn. 50) is linear, it will benefit from expansion of \G) in the eigenfunctions \\j) of 
the operator C, where 

t\X i )=X j \X i ) (62) 

for eigenvalues Xj. In Appendix G we show that the commutation relation (Eqn. 61) and the steady state solution 
(x\G) = e 9 ( x_1 ) , for which C\G) — 0, completely define the eigenfunctions and eigenvalues of C; we summarize the 
results of Appendix G here. The eigenvalues are shown to be nonnegative integers j, 

Xj=j€ {0,1,2,...}, (63) 

as in Eqn. 42. The eigenvalue equation now reads 

qj)=b + b-\j)=j\j), (64) 

which is consistent with interpretation of b + b~ as a number operator for the eigenstates \j). The eigenfunctions are 
shown to be (in x space) 

(x\j) = (x-iye g ( x - 1 \ (65) 



as in Eqn. 43 with <f>j(x) = (x\j). The conjugate eigenfunctions are shown to be 

-g(x-l) 

m - (^dTTI- (66) 

The operators b + and b~ are shown to raise and lower the cigcnstates \j), respectively, as a + and a - do the \n) states, 
i.e. 

b + \j) = 1.7 + 1), (67) 

b~\j) = j|j — 1>, (68) 

(j\b + = 0'-l|. (69) 

(j\b- = (j + l)0' + l|. (70) 

Now employing the expansion of \G) in the cigcnstates \j), 

\G(t)) = J2 G i®W> ( 71 ) 

3 

the equation of motion (Eqn. 50) gives a trivial equation for the expansion coefficients 

Gj = -jGj, (72) 

which is solved by 

G 3 {t)=A je -^ (73) 

for some constants Aj, as in Eqn. 40. The inverse of Eqn. 71 (obtained by projecting (j\ and using (j\f) = Sjj 1 ) is 

U\G(t))=Gj(t). (74) 

The probability distribution is retrieved by inverse transform (Eqn. 49), 

Pn (t) = (n\G(t)) = (n^GMJ) = X^e-^ntf), (75) 

where the coefficients ^4j are computed from the initial distribution p n (0): 

4* = Gi(0) = (j|G(0)) = J2p n (0)(j\n). (76) 

n 

Eqns. 75 and 76 give the full time-dependent solution to the birth-death process as expansions in the eigenmodes 
(n\j) and conjugate eigenmodes (j\n). Because we have decomposed the problem using the eigenbasis, or spectrum, 
of the underlying operator, we will refer to this as the spectral solution. 

The quantities (n\j) and (j\n) are overlaps between the protein number basis \n) and the eigenbasis \j). Both (n\j) 
and (j\n) are readily computed by contour integration or more efficiently by recursive updating; these techniques are 
presented in Appendix H. Notable special cases are 

H0> = e-° 9 - (77) 

n! 

(0|n) = 1, (78) 

which confirm that Eqn. 75 describes a Poisson distribution in steady state (j = 0). 

C. Fokker-Planck approximation 

The previous sections discuss several methods for calculating the steady state and time dependent solutions of the 
master equation for the birth-death process. For many larger systems with regulation, full solution of the master 
equation is not possible. In some of these cases one can make progress by deriving an approximate equation which is 



10 

valid when protein numbers are large. In the limit of large protein numbers the master equation can be expanded to 
second order to yield the Fokker-Planck equation. In this section we derive the Fokker-Planck equation for a general 
one-dimensional master equation with arbitrary production and degradation rates; the method is easily generalizablc 
to models with regulation. 

For arbitrary production rate g n and degradation rate r n , the master equation reads 

OtPn = - (g n + r n )Pn + g n -lPn-l + r n +lPn+l 5 (79) 

setting g n = g and r n = n recovers the simple birth-death process with time rescaled by the degradation rate, as in 
Eqn. 9. 

1. Large protein number 

The Fokker-Plank equation is derived under the assumption that the typical protein number is large (n» 1), such 
that n can be approximated as a continuous variable, and a change of 1 protein can be treated as a small change. 
We will use parentheses when treating n as continuous and a subscript when treating n as discrete. Under this 
approximation, the function f(n ± 1), where f(n) <G {g{n)p(n),r(n)p(n)}, can be expanded to second order as 

f(n ± 1) = f(n) ± d n f(n) + l -d 2 J(n). (80) 

Inserting the results of the expansion into Eqn. 79, we obtain 

d t p(n) = -d n [v(n)p{n)\ + ^d 2 n [D(n)p(n)\ , (81) 

where v(n) = g(n) — r(n), an effective drift velocity, recovers the right-hand side of the deterministic equation, and 
D{n) = g(n) + r(n), an effective diffusion constant, sets the scale of the fluctuations in protein number. The drift 
term plays the role of an effective force. The diffusion term plays the role of an effective temperature: the larger it 
is, the larger the excursions a single trajectory (of particle number versus time) takes from the mean. Eqn. 81 is the 
Fokkcr - Planck equation. 

The steady state solution of Eqn. 81 is easily obtained by noticing that the Fokker-Planck equation is a continuity 
equation of the form 

dtp = -d n j, (82) 

where j(ri) = v(n)p(n) — \d n [D(n)p(n)} is the current of the probability. In steady state (d t p = 0) the current is 
constant, and since p(n) — > and d n p(n) — ^ as n — > oo (typically more quickly than v(n) or D{n) diverges), the 
current vanishes for all n. The steady state distribution is then found by direct integration, i.e. 



< \ N 
P{n) = D{n) 6XP 



dn' v{n,) 



D(n') 



(83) 



where N is a normalization constant ensuring J„ dnp(n) = 1. 

In the simple birth-death process, for which v(n) = g — n and D(n) = g + n, Eqn. 83 evaluates to 

p(n) = - ( 1 + - V ^ 2n - (84) 

In Appendix D, we show that the exact steady state (the Poisson distribution, Eqn. 16) and Eqn. 84 have the same 
asymptotic limit for large protein number: a Gaussian distribution with mean and variance equal to g, the ratio of 
production to degradation rates. 

2. Small noise 

In addition to the approximation that the typical protein number is large, one may make the further approximation 
that the noise is small. This is often referred to as the "linear noise approximation" [17, 37, 38] or "small noise 
approximation" [32, 33]. Specifically, one assumes that the fluctuations rj in n are small around the mean n, i.e. 



n + r], (85) 



11 

where n = (n) = ^2 n np n (for brevity we use bar notation in this and several subsequent sections) . Since dn/dr/ = 1 
we have p(n) = p(ji), and the master equation (Eqn. 81) becomes an equation in rj: 

dtp(v) = -d v [v{n + v)p(v)} + ^ d v l D ( n + rfPiv)] , (86) 

We use the approximation that rj is small to expand the drift and diffusion terms to first nonzero order in 77: 

v(n + r/) = v(n) 4- rjv (n) + ■ ■ ■ « rjv (n) (87) 

D(n + r]) = D(n)+77L>'(n)H « D(n) 

where prime denotes differentiation with respect to n, and the last step in Eqn. 87 recalls the fact that v(n) 
g(n) — r(n) = in steady state, as given by the kinetic rate equation (e.g. Eqn. 8). Eqn. 86 is now 



d t p(v) = -v'(n)d v [m>(ri)\ + -D(n)d 2 \p{rj)\ . 



As with Eqn. 83, the steady state of Eqn. 89 is found by direct integration: 



p{rj) 



where 



N 


\ 2V ' {fi) f eft/V 

[D(n) Jo V \ 


1 


n/ _s cxp 
D(n) 


V27TCT 2 




2 D(n) 
a = — ——. — r 





V/2^ 2 



-2v'(n) 



(89) 



(90) 



(91) 



(note that v'(n) is negative for stable fixed points). Eqn. 90 is a Gaussian distribution, meaning that under the linear 
noise approximation, the steady state probability distribution is a Gaussian centered at the exact mean n and with 
width determined by mean birth and death rates according to Eqn. 91. We note that Eqn. 90 can be equivalcntly 
derived by expanding the integrand in Eqn. 83 about its maximum (or "saddle point") n. 

The linear noise approximation is stricter than the large-protein number approximation made in the previous section 
(Sec. IIC 1). While the previous approximation makes no assumption about the form of the probability distribution, 
the linear noise approximation assumes that the distribution is unimodal and sharply peaked around its mean. In 
practice, n and a 2 are obtained by finding the steady state(s) (i.e. the stable fixed point(s)) of the corresponding 
deterministic rate equation. However, it is easily possible (for processes more complicated than the simple birth-death 
process) for the deterministic equation to have more than one stable fixed point (see Fig. 2). Although one may make 
a Gaussian expansion around each fixed point in turn, the linear noise approximation does not describe how these 
Gaussians might be weighted in a multi-modal distribution. In these cases it is most accurate to use (if solvable) 
cither the large-protein number approximation (Eqn. 83) or the original master equation (Eqn. 79). 

In the simple birth-death process, for which v(n) = g — n, D(n) = g + n, and n = g, Eqn. 91 gives a 2 = 
2g/[— 2(— 1)] = g, and therefore Eqn. 90 reproduces the asymptotic behavior derived in Appendix D. 

D. Langevin approximation 

We now consider a second stochastic approximation: the Langevin equation. The advantage of the Langevin 
approach is that a large amount can be learned about the process without finding the full distribution, but instead 
by considering the correlation functions of the concentration, which arc readily computed. Starting from the Fokker- 
Planck equation we can calculate the equation for the mean of the distribution (n) = ^2 n= Qnp n . We arrive at the 
kinetic rate equation given in Eqn. 2. 

We can think about the change in concentration in time as a trajectory in n space. Each realization of a birth-death 
process, will be described by a certain n(t), and averaging over many experiments, we obtain (n(t)), as given by Eqn. 
2. If we consider the change of concentrations on time scales longer than the characteristic timescales of the particular 
reactions, we can assume that there are many birth and death processes in each interval and that the fluctuations of 
each realization of n(t) around the mean are Gaussian distributed: 



P\rj\ 



■Jdnfdt^il 



(92) 
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We can include fluctuations around the the mean values by considering an additional noise term f](t), such that the 
equation for the change of the concentration of proteins n becomes: 

U/Tl 

— = v(n) + r](t) = g-n + r)(t). (93) 

at 

We require of the noise that: 

(V(t)) = 0, (94) 

(vtf)v(t)) = S{t-t')D{n) = 5(t-t')(g + n) (95) 

where D(n) is in general the diffusion term in the Fokker-Planck equation. In Appendix I we show the equivalence 
of the Fokker Planck and Langevin descriptions. 

In the general case when there are many types of proteins in the system, we can define a time dependent correlation 
function 

Cij(t) = (SNiQtySNjit)), (96) 

where the average implies a time average and SNi(t) = Ni(t) — (6Ni(i)) are the variances of each type of protein. In 
one dimension Sn 2 is the variance a 2 (in subsequent sections we use Sn 2 and a 2 interchangeably). We note that in 
steady state the time average can be replaced by an ensemble average. For the case of the single species birth-death 
process we have already computed the means and variances in Sec. II Bl. In case of the birth-death process, we 
calculate the time dependent correlation function to obey the equation: 

In steady state the solution must reduce to the previously calculated variance (Appendix C): Sn 2 = 
((n — (n))((n — (n)))) = (n 2 ) — (n) 2 = g. Therefore the solution of Eqn. 97 for a simple birth-death process is 

C(t) = ge-K (98) 

We can also consider the correlation functions in Fourier space 

C,y(w) = f dte-^ l C l0 {t). (99) 

Jo 

Often it is easier to calculate the Fourier transform of the correlation function directly from the Fourier transform of 
the Langevin equation, and then invert the transform back. In case of the simple birth-death process we can vary the 
Langevin equation around its mean values n(t) — (n) + Sn(t) (that is simply linearize the equation around its mean) 
and consider the resulting equation in Fourier space to obtain 

— iu>6h(ui) = —5n{uj) + fj(ui), (100) 

where 5h{ui) = J °° dte~ luJt n{t) and 

/>oo 

fj(u) = / dte- iwt r](t) (101) 

Jo 

(jKw)tKu/)) = 2ir5(Lu-io')(g+(n)). (102) 



The correlation function is then: 



(6n*(uj)Sn(u')) = 2tt5(lj - u') - 't — = 2ir^—. (103) 

U1 Z + 1 LO Z + 1 

Inverting the Fourier transform reproduces the real time correlation function. 

C(t) = r due iu}t ^— = 2Tne^ = ge~*, (104) 

Jo lu 2 + 1 2i v ' 

where we reproduce the result of Eqn. 99. The real part of the auto correlation function in Fourier space is also called 
the power spectrum 

/>oo 

Af(u) = (Sn*(w)Sn(u/))= dte- lU)t C ll {t), (105) 

Jo 

because it tells us which frequency modes contribute most to the form of the noise. We note that the power spectrum 
N(lo) may be written as integral of the correlation function in real time. 
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E. Comparison of the descriptions 

For the birth-death process, the steady state of the kinetic rate equation agrees with the mean of the steady state 
probability distribution (under both the exact and approximate stochastic descriptions): (n) = g. In the limit of large 
protein number, the Poisson distribution is well approximated by a Gaussian, and the Fokker-Planck (and therefore 
the Langevin) approximation is a good description. To investigate the validity of the Fokker-Planck approximation 
at small and intermediate protein number, in Fig. 1 we compare for a range of protein numbers the probability 
distribution obtained directly from the master equation (i.e. the Poisson distribution, Eqn. 16) and from the Fokker- 
Planck approximation (Eqn. 83). We quantify the disagreement between two distributions using the Kullback-Leiblcr 
divergence, 

£>XL = $>„log^, (106) 

where p n corresponds to the master equation and p n corresponds to the Fokker-Planck approximation. The Kullback- 
Leibler divergence is not symmetric and is therefore appropriate for a comparison between a "true" distribution (p„ 
here) and its approximation (p n here). As the mean protein number (n) = g increases, the accuracy of the approxi- 
mation increases, and the divergence Dkl decreases. We plot explicitly three sample pairs of exact and approximate 
distributions, at small, intermediate, and large protein numbers. As expected the Fokker-Planck distribution deviates 
from the Poisson distribution at small protein number, and agrees well at large protein number. 

III. AUTOREGULATION 

We now begin to turn our attention from the simple one-dimensional birth-death process to more realistic models 
gene regulation (see [43, 44] for a discussion of regulation functions of gene expression). We start with the description 
of auto-regulation of one gene, in which its protein production rate is an arbitrary function of its own protein number. 

A. Deterministic model 

As before the mean dynamics are captured by the kinetic rate equation. The kinetic rate equation for a birth death 
process in which the production rate is an arbitrary function g{n) of mean protein number n is 

dn , „ ... _. 

— = g(n) - n. (107) 

The problem becomes potentially much harder, since the auto-regulation function g{n) can be nonlinear. 

The form of the auto-regulation function depends specifically on the molecular model of regulation which is being 
considered. For example, the Hill function, 

g_K h +g+n h 
9{n) = K* + n* ' (108) 

can be derived by considering a gene with two production rates g+ and <?_ (corresponding to the states in which a 
protein is bound and unbound to the DNA, respectively) for which the rate of switching to the bound state depends 
on the protein number (see Sec. IV and Appendix J). The derivation assumes that the rates of binding and unbinding 
are faster than the protein degradation rate (which sets the time-scale for changes in protein number), such that 
equilibrium is reached, with equilibrium binding constant K. 

The parameter h describes the cooperativity of protein binding. When g + > g_, the case h > corresponds to 
activation and the case h < corresponds to repression; the special case h = reproduces the simple birth-death 
process. For \h\ > 2, Eqn. 107 has two stable fixed points for certain parameter regimes; for \h\ > 2 fixed points must 
be found numerically. 

Although the Hill equation is often used in models of gene regulation, we note that other functional forms are 
derivable from other biochemical processes. Many of the following results are valid for arbitrary g(n). 

B. The master equation 

The full stochastic description corresponding to the deterministic Eqn. 107 is given by the birth-death master 
equation in which the production rate is described by the arbitrary auto-regulation function g n (recall that we replace 
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FIG. 1: Comparison of the distributions obtained from the master equation (Eqn. 16) and from the Fokker-Planck approx- 
imation (Eqn. 83), for the simple birth-death process. The Kullback-Leibler divergence Dkl (Eqn. 106) between the two 
distributions is plotted as a function of mean protein number (n) = g. Insets show distributions obtained from the master 
equation (solid) and from the Fokker-Planck equation (dashed) for g — 0.1 (star), g — 0.9 (circle), and g — 20 (cross), where 
the symbols correspond to points on the Dkl curve. 



parentheses with subscript when treating n as discrete): 

Vn = -gnPn ~ np n + g n -ip n -i + (n + l)p n +l, 

We may easily generalize the solution in Sec. II B 1 to find the steady state probability distribution, 



Pr, 



P() 



n-1 

n 



(109) 



(110) 



with po set by normalization. Except for special cases of the regulation function we cannot find a closed form solution 
for the distribution, but the product is easily evaluated. 

The results in Sec. II for the simple birth-death process can be generalized in the case of auto-regulation. We will 
pay particular attention to the generalization of the eigenfunction expansion in operator notation for arbitrary g n , 
when discussing two genes in Sec. V. 
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C. Bistability and noise 

Auto-regulation can affect the statistics of the steady state distribution. Even before specifying the form of g(n), 
we may obtain a general statistical result in the limit of large protein number by using the linear noise approximation 
(Sec. IIC 2). Eqn. 91 describes the variance a 2 of fluctuations around the steady state mean n. Here we compute the 
ratio a 2 jn for auto-regulation in order to compare with the Poisson distribution, for which a 2 /n = 1. 

For auto-regulation the diffusion term becomes D{fi) = g(n)+h = 2n, where we have used the fact that = g(n) — n 
at steady state. The derivative of the drift term evaluated at the mean is v'(n) = g'(n) — 1, where the prime denotes 
differentiation with respect to n. Thus Eqn. 91 becomes 

fj2 - 1 

~n ~ l-ff'(n)' (111) 

This expression shows that self-activation (g'(fi) > 0) leads to super- Poissonian noise (a 2 jn > 1), and self- repression 
(g'(n) < 0) leads to sub-Poissonian noise [a 2 jn < 1); when there is no regulation we have g'(n) = 0, and we recover 
the Poisson result, a 2 jn = 1. Fig. 3 demonstrates this behavior by computing the exact distribution (Eqn. 110) for 
Hill function auto-regulation (Eqn. 108) with various values of the cooperativity parameter h. 

Eqn. Ill diverges when g'(n) = 1. This corresponds to a bifurcation point for the self-activating gene, where there 
are two possible solutions to the steady state equation. In this case the ratio of the variance to the mean is not very 
informative about the distribution because the distribution is bimodal. The self-repressing gene, in contrast, cannot 
be bistable, and the ratio of the variance to the mean is always a good description of noise. In Fig. 2 we demonstrate 
that the exact distribution becomes bimodal when the deterministic equation crosses a bifurcation point. 

Eqn. Ill can be be equivalently obtained by computing the power spectrum from the linearized Langevin equation. 
Linearizing the Langevin equation and going to Fourier space we obtain 

— iu)5h(ijj) = g' (h)h(uj) — n(ui) + fj(u>). (112) 

The power spectrum becomes 

2n 
Co 2 + [1 - g'(n)} 2 v ' 

The steady state auto-correlation function is therefore 

Sn 2 = n (114) 

\l-g'(n)\ 



and the ratio of the variance to the mean is 

Sn 2 



n |1 — g'(n) 



(115) 



We point the reader to the paper of Warren, Tanase-Nicola and ten Wolde [35] for a pedagogical discussion of power 
spectra in gene expression. 

IV. BURSTS 

The previous section assumes that the processes of transcriptional and translational regulation can be described 
by one deterministic regulation function. Recent experiments have observed a feature of gene regulation that is not 
captured by such a model: protein production often occurs in bursts [4, 45, 46]. Bursts can result from the gene 
transitioning among two or more states with different rates of protein production, for example when a transcription 
factor is bound or unbound, or, in higher eukaryotes, when different production states are introduced by chromatin 
remodeling [45] . Bursts can also arise during translation, in which a single mRNA transcript produces many copies of 
the protein in a short amount of time. Here we present simple models of the first type, which we call transcriptional 
bursts, and the second type, which we call translational bursts. We then discuss a model which combines bursting 
with auto-regulation. We refer the reader to the literature for discussion of other burst models including those which 
capture both transcriptional and translational bursts [25, 47]. 
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FIG. 2: The steady state solution(s) for a self-activating gene as a function of the basal rate <?_. The figure shows regimes 
which are mono-stable (high and low gr_) and bistable (intermediate <;_). The probability distributions below (solid black lines) 
correspond consecutively to the points above marked by stars. In the mono-stable regime the distributions are compared to 
Poisson distributions (solid gray lines). The auto-regulation function is a Hill function (Eqn. 108) with h — 4, g+ — 10, and 
K = 6. 

A. Transcriptional bursts 

In this section we consider a gene that can exist in multiple states, each with its own protein production rate. 
We first specialize to the case of two states, which is commonly used as a model for bursting [4, 34]. We solve this 
case conventionally using a generating function. The conventional solution is limited to two states; for an arbitrary 
number of states we present the spectral solution using operator methods. Finally we present a solution to a simple 
model with both bursting and auto-regulation. 



The master equation 



We begin by considering the general process in which a gene has Z different states. In each state z the protein 
production is described by a simple birth-death process (Sec. II) with production rate g z . Transitions between the 
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FIG. 3: The ratio of the variance Sn 2 to the mean (n) as a function of the cooperativity coefficient h, for a gene undergoing 
Hill function (Eqn. 108) autoregulation. For h < the gene represses its own production; for h > the gene activates its own 
production. The case h — recovers the simple birth-death process result Sn 2 /(n) — 1. Points are plotted for integer values of 
h; the line shows Eqn. 108 for continuous h. Insets show probability distributions (black solid lines) from the master equation, 
compared with Poisson distributions with the same mean (gray solid lines), at h — — 1 (star), h — (circle) and h = 1 (cross), 
where symbols correspond to points on the 8n 2 /(n) curve. Parameters are <?_ = 2, g+ — 20, and K — 4, which restrict solutions 
to the mono-stable regime, yielding unimodal distributions. 



states occur at rates given by the transition matrix fl ZZ ', such that 

z 



E n « 



0, 



(116) 



which follows from conservation of probability. The rates are constant here; in Sec. IV A 4 we extend the model such 
that the rates depend on protein number. 

The system is described by a Z-state probability vector with elements given by p„- The particular states evolve 
according to the coupled master equation 



V'n 



-9zPn - np z n + g z p z n _ x + (n + l)p z n+1 + ^ ^zz'Pn ■ 



(117) 



where all but the last term describes the birth-death process for state z, and the last term describes transitions among 
states. The probabilities of being in state z regardless of the number of proteins are given by 7r z = ^2 n Pn', summing 
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the master equation over n in steady state gives 

5^n„/7r z /=0, (118) 

z' 

and normalization requires 

X>* = i- ( 11Q ) 

z 

2. The two-state gene 

We first consider the case in which the gene has only two states (Z = 2), an active state (z = +) and an inactive 

state [z = — ): 

A = (*) ■ (120) 

The transition matrix takes the form 

n ~' = {7: X) • (121) 

where lo + and w_ are the transition rates to the active and to the inactive states, respectively. 

As in Sec. II B 2, we define the generating function G±(x) = J2 n Pn xn an d sum the master equation (Eqn. 117) 
over x n , yielding at steady state 

= -y(dy-g±)G ± ±uj + G-Tu-G + , (122) 

where y = x — 1. Writing G±(y) = e 9±y H±(y) yields a simpler equation for H±(y); writing out the + and — cases 
explicitly, 

= -ye yA d y H + +L} + H--LO-e yA H+, (123) 

= -ye- yA d y H-+oj-H + -Lj + e- yA H-, (124) 

where A = g + — g_. Eqns. 123-124 are a coupled pair of first-order ordinary differential equations; combining them 
yields a second-order differential equation. Specifically, solving Eqn. 123 (124) for H- (H + ) and substituting it into 
Eqn. 124 (123), one obtains 

= uc%H± +(/3- u)d u H ± - aH±, (125) 

where u = +yA, a = W:p, and ft = w + + w_ + 1. Eqn. 125 is the canonical equation for the confluent hypcrgcomctric 
function. Only one of its two solutions satisfies p n — > as n — > oo; it is H±(x) = N±Q[a, ft; u], where 

^r a i ^r(^ + a ) r(/3) u e , 10C . 

$ K^] = EViW)I (126) 

is the confluent hypergeometric function of the first kind, and the scalars N± are found by normalization: N± — 
G±(l) = J2 n Pn( 1 ) n = v ± = lu±/(lj+ +LJ-) (the last step uses Eqns. 118-119). Thus, 

G±{x) = "* e »^ x -V*[LJ T ,u> + + w_ + 1; T(g+ - 9-)(x - 1)]. (127) 

As before, the probability distribution can be recovered from G(x) = X]± G±(x) by differentiation (Eqn. 21). 
As shown in Appendix K, in the limit <?_ — 0, Eqn. 127 reduces to the slightly simpler expression 

G(x)=*[w + ,w++w-;g + {x-l)], (128) 

for which the probability distribution has the explicit form 

g+T(n + uj + ) r(w + +w_) . n ,„ . 

n\ 1 \u>+) l(n + w + +w_) 

in agreement with the results of Iyer-Biswas et al. [34] and Raj et al. [4] . 
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3. The spectral solution 

It is clear that the method presented in the previous section results in a Zth order ordinary differential equation, 
and therefore its utility is generally limited to Z = 2 states. In contrast, here we show that the spectral solution, in 
which we expand in eigcnfunctions of the birth-death operator, can be used for an arbitrary number of states. 

We begin as in Sec. IIB3 by defining \G Z ) = J2 n Pn\ n ) an d writing the master equation in terms of raising and 
lowering operators, 

\G z ) = -b+b-\G z )+Y,n zz '\G z ,), (130) 

z' 

where now b~ = a~ — g z is z-dcpcndent. We expand the generating function in the z-dcpcndent eigcnfunctions of the 
birth-death process as 

|G 2 ) = ^G||i,), (131) 

3 

where the eigenvalue relation reads 

b + b-\ Jz )=j\ ]z ). (132) 

This gives the following equation for the expansion coefficients: 

G ] = -3G] + £ n zz , J2 Of ( Jz \j' z> ). (133) 

z' j' 

Here the transition matrix £l zz > couples the otherwise independent birth-death processes, and 

(jz\j Z ')= { ~ AzZ ' y ,V 0(j-f), (134) 

u -fr- 

where A zz > — g z — g z i, and the convention 8(0) = 1 is used for the Heaviside function. Eqn. 134 is derived by 
evaluating the inner product using the x space representations of the eigcnfunctions, 

(x\j z ) = (x-l) i e<"<- x - 1 \ (135) 

g-9zO-l) 

Ux\x) = {x ~ 1)J+1 , (136) 

and Cauchy's theorem, as in Appendix E. 

The Heaviside function makes Eqn. 133 lower-triangular. This is explicitly clear, for example, in steady state, 

jo] -zZ^' G 1 - E a " E G i ( "rf r'iV ' (137) 

where one observes that the jth component is a function only of components j' < j. The lower-triangular structure 
allows G z z to be computed iteratively, which makes the spectral solution numerically efficient. The lower-triangular 
structure is a consequence of rotating to the eigenbasis of the birth-death operator; this structure was not present in 
the original master equation. Other spectral decompositions are possible by exploiting other eigenbases; for example 
by expanding in a single eigenbasis parameterized by a constant (not z-dependent) production rate, one obtains an 
equation that is sub-diagonal in j [29] . 

The probability distribution is recovered via inverse transform, 

ti = J2 G 'MJ*)> ( 138 ) 

j 

as in Sec. II B 3. The spectral solution is valid for any number of states Z, while direct solution of the differential equa- 
tion is generally limited to Z — 2. As expected, for Z = 2 the spectral solution recovers the confluent hypcrgeometric 
function solution (Eqn. 127), as shown in Appendix K. 

Fig. 4 demonstrates the spectral solution for Z = 3 states when the transition rate among states is (a) equal to, and 
(b) much less than, the degradation rate. In (a) the bursts give rise to a long tail in the distribution (compared to 
the Poisson distribution) ; in (b) the system spends long times in each production state before transitioning, resulting 
in a trimodal distribution. 
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FIG. 4: Steady state probability distribution for transcriptional bursting with Z = 3 states, when the transition rate among 
states is (a) equal to the degradation rate (oi = 1), and (b) much less than the degradation rate (uj — 0.01). Production rates 
are g z = (0, 3, 15) and the rate of transition u> between any two states is equal (i.e. all off-diagonal elements of £l zz t are u> and 
all three diagonal elements are — 2uj). 



4- Bursting and auto-regulation 

Here we consider the simplest model that combines bursting with auto-regulation [27]. The solution follows that of 
Sec. IV A 2. The gene has two states, an active state in which the regulatory protein is bound (p^) and an inactive 
state in which the regulatory protein is unbound (p~). The regulation is incorporated by making the rate at which 
the gene switches to the active state proportional to the protein number: 



n xr/ = 



— cj_ uj + n 



(139) 



In this case the analog of Eqn. 122 for the generating function is 

= -(a; - l)(d x - g±)G± ± u+xd x G- T U-G+, (140) 

As in Sec. IV A 2, the + and — equations can be combined to give a second-order differential equation, e.g. for G_(x), 

= d x ;G- + C 1 (x)d x G-+Co(x)G-, (141) 

with 

(142) 

(143) 

whose solution G-(x) = iV_e 5,+ ( 2:_1 '$[a,/3;u] is proportional to the confluent hypergeometric function of the first 
kind (Eqn. 125) with 



Ci(a;) = 
C 2 (x) = 


9- 


-9+x- g. 

(1 + 


(1 + 
- (.9+ + w_ 
d + )x — 1 


+ 1) 


-1 



a = 1 

(3 = 1 



l + w+ 

U>- 



w+ff- 



9- - (l+^+).9+ 



1+UJ+ (l+w+) 2 ' 
[. 9+ (l+^ + )-ff_][(l+^ + )x-l] 

(l+^+) 2 



(144) 
(145) 

(146) 



21 

Again 7V_ is set by normalization; G+(x) is computed from G-(x) using Eqn. 140 (bottom signs). 

In Appendix J we consider an extension of this model, in which the transition rate to the active state is proportional 
to n , where h is the cooperativity of the binding. We show that in the limit of fast transitions this model yields an 
effective auto-regulation function of the Hill form (Eqn. 108). Full solution of this model using the method presented 
here for h = and h = 1 is of limited utility for h > 1, because the high-order differential equation is not solvable. 



B. Translational bursts 



We now turn our attention to translational bursts, in which a single mRNA transcript produces many proteins 
in a short amount of time. In this case the model includes just one protein production rate, but each production 
event results in an instantaneous burst of multiple proteins instead of one. This approach does not explicitly consider 
mRNAs; instead it models the translational step as an effective burst of proteins. 

With production occurring in bursts of size N, the birth-death master equation becomes 



P„ 



-gpn 



nPn + gPn-N + (n + l)Pn+l 



(147) 



This equation is easily solved by spectral decomposition. In terms of raising and lowering operators (Sec. II B 3), Eqn. 
147 reads 



\G) 



, A r 



-g-a + a + g(a + ) +a \G) 

-b+b- + r)\G), 



(148) 
(149) 



>iV 



where b + = a + — 1 and b = a — g as before, and T = (a + ) — a + . Just as in Sec. IIB3, Eqn. 149 benefits from a 
spectral expansion in the eigenfunctions \j) of the simple birth-death process, i.e. 



\G) = J2 G M' 



(150) 



where 

b + b-\j)=j\j). 
Eqn. 149 then becomes an equation for the expansion coefficients: 

Gj^-jGj + YsdimGj 



(151) 
(152) 



The simplification 



(j\m 


= 01 


(V + i^-fV + i) 


\f) 




= 01 


.£=0 v / 




= 01 


£=2 ^ ' 


r) 



N 



= (^-i)^'-i+E 



1=2 



N 
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N 



G, = -jG, + (N-l)G^ 1+ J2 



N 



e=2 



(153) 
(154) 

(155) 

(156) 



(157) 
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The equation is lower-triangular, which makes a recursive solution simple. Note that setting N = 1 recovers Eqn. 72 
for the simple birth-death process. The probability distribution is retrieved by inverse transform, 



P, 



= ]TG>|j). (158) 



with (n\j) computed as in Appendix H. 

Translational bursts increase noise relative to the simple birth-death process. For example, we can easily compute 
the variance and mean of the steady state distribution for translational bursts (see Appendix C); their ratio is 

a 2 1 + N , N 

¥ = — < 159 > 

which for N > 1 is (potentially much) greater than the birth-death result a 1 jn = 1. 
We note that within the Langevin approach we can model translational bursts as: 

dn , „ , 

— = Ng-n + r](t) (160) 

at 

(v(t)v(t')) = 5{t-t')(Ng + n). (161) 

We obtain a power spectrum of the form: 

AT=^±i (162) 

u> z + 1 

which again is a signature of a more noisy process. 

V. TWO GENES 

At this point we have discussed models of a single gene, including the simple birth-death process, auto-regulation, 
and bursts. Here we extend the discussion to include multiple genes with regulation. We focus on the case of two 
genes, and we describe how the method may be extended to larger regulation networks. We refer the reader to the 
literature for more detailed discussion of the application of the spectral method to cascades [28] and to systems with 
both regulation and bursts [29]. 

Here we consider two genes, the first with protein number n, and the second with protein number m. The first gene 
regulates its own expression via arbitrary function g ni and it regulates the expression of the second gene via arbitrary 
function q n . This system can be thought of as the simplest possible regulatory network; it is also the minimum system 
necessary for completely describing a regulatory cascade of arbitrary length [28] . Because the model is general, the 
two degrees of freedom could correspond to two species of proteins, or to a protein and an mRNA species. 

The master equation for this system is 

Prim QnPnm ^iPnm i Qn—lPn—l,m i \7l T" -L/Pn+l,m 

+P [-QnPnm - mp nm + q n p n ,m-l + (m + l)p n ,m+l] , (163) 

where time is rescaled by the degradation rate of the first gene and p is the ratio of the degradation rate of the first 
gene to that of the second. 

A. Deterministic model 

As always, the kinetic rate equations can be obtained by deriving the mean dynamics from the master equation. 
Summing the master equation over both indices against either n or m (and performing index shifts as in Appendix 
B) gives 

dn _ 

-jr = 9\n) - n, (164) 

dt 

— — = q(n) — fh, (165) 

dt 
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respectively. Here we have assumed that ^ n g n p n w g(n) and J^ n g n p„ w q{n) by taking n = n + 77 and keeping the 
first terms of the Taylor expansions g(n) — g(n) + r]g'(n) + . . . and q(n) = q{n) + r)q'{n) + . . . , respectively. 

Solving for the steady state of the deterministic system requires finding the roots of two coupled nonlinear equations, 
which can only be done numerically Solving the deterministic system is essential for finding the power spectrum and 
calculating the correlation functions in the Langevin approach. However, as we show in the next section, the full 
stochastic description can be solved quite efficiently using the spectral method, i.e. by expanding in the cigenfunctions 
of the simple birth-death process. 

B. The spectral method 

Although it is usually not possible to find an exact solution at either the master and Fokker-Planck equation levels, 
here we show that efficient numerical computation of the full steady state distribution is possible by expanding in the 
cigenfunctions of the simple birth-death process (Sec. II B 3). This approach is general and does not depend on the 
algebraic forms of the regulation functions. 

As in Sec. IIB3 we define the generating function \G) — ^2 nm Pnm\n,m) , which makes the master equation 

\G) = -£\G), (166) 

where the full linear operator 

£ = b+b-(n)+pb+ n b m (n) (167) 

is written in terms of raising and lowering operators defined analogously to Eqns. 59-60: 

b+ = at - 1, (168) 

b+ = a+ - 1, (169) 

b~(n) = a~ -g{n), (170) 

&"(n) = a m -q(n). (171) 

The difference here is that both lowering operators are n-dependent because of the regulation functions. This suggests 
that we separate the full operator as C = £0 + Ci(n), where we have defined a constant part, 

£ = 6+6" + pS+6" , (172) 

with b~ = a~ — g and 6~ = a~ — q, and an n-dependent part 

A(n) = 6 + f(n) + pS+A(n), (173) 

with r(n) = g — g(n) and A(n) = q — q(n). The operator Cq describes two independent birth-death processes with 
constant production rates g and q respectively, and the operator C\ captures how the regulated processes differ from 
the constant rate processes. The constants g and q act as gauges: they can be chosen freely, and the final solution is 
independent of this choice. 

We have already solved for the eigenfunctions |j, k) of the constant rate birth-death processes in Sec. II B 3. The 
spectral method exploits this fact by expanding the solution of the full problem in these eigenfunctions: 

\G)=Y,G lk \j,k). (174) 

Eqn. 166 then gives an algebraic relation among the expansion coefficients: 

Gjk = -(.?' + pk)G jk - Y^ T j-ij'Gj<k -p^Ajj.Gj^k-!, (175) 

i' i' 



where 



r,,< = {mn)\j') = Y,{J\n){9-9n){n\j'), (176) 

n 

A ]r = (j\A(n)\f) =^(j\n)(q- q n )(n\f). (177) 
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FIG. 5: Demonstration of the spectral method with two genes. Regulation is of Hill form, q n — (q-K h + q + n h ) / (K h + n h ), 
with q- — 1, q+ — 12, K — 4, and h — 4. In (a) the first gene undergoes a simple birth-death process with g„ = 7; in (b) the 
first gene regulates itself with g„ — q n . Degradation rates are equal: p — 1. 



Eqn. 175 enjoys the property that it is sub-diagonal in the second gene's eigenvalue k. This feature conies from the 
fact that both regulation functions are dependent on the first gene's protein number n (and not m). The sub-diagonal 
structure allows the matrix Gjk to be computed iteratively by columns in k, which makes the method numerically 
efficient. The initial condition 



G j0 = (j, k = 0|G) = 5>nm0»(0|m) = Y,Pn(j\n) 



(178) 



requires the marginal probability distribution of the first gene p n , which can be found recursively (Eqn. 110) as 
described in Sec. Ill B on one-dimensional auto regulation. The joint probability distribution is retrieved by inverse 
transform: 



Pnm = J~]Gjk(n\j)(m\k). 



jk 



(179) 



To demonstrate the capability of the spectral method, Fig. 5 plots the steady state joint probability distribution 
for a particular choice of regulation function, for cases (a) without and (b) with auto-regulation. To demonstrate 
the accuracy of the spectral method, Fig. 6 compares steady state marginal distributions obtained by the spectral 
method with those obtained both by direct iterative solution of Eqn. 163 and by standard Gillespie simulation [48]. 
The spectral method attains excellent agreement with both methods in orders of magnitude less runtime than either 
other method; full details on the efficiency and accuracy of the spectral method are available in [28]. 

The spectral decomposition presented here constitutes an expansion in only one of many possible eigenbases: that 
with constant production rates. Other spectral expansions are possible, as discussed in detail in [29]; we do not repeat 
the discussion here. 



VI. SUMMARY 



In this pedagogical chapter, we have attempted to introduce the reader to the phenomenon of intrinsic noise and 
to the simplest possible mathematical tools for modeling. We have avoided detailed discussion of a vast literature 
in Monte Carlo simulation of such phenomena in favor of simple calculations which can be performed analytically 
or which are useful in numerical solutions of the master equation itself (rather than random generation of sample 
trajectories). 

The discussion also hopefully illustrates how different problem settings motivate different analytic methods, many of 
which have close parallels in the physics canon. These include approximation of differences with derivatives, yielding 
a description very similar to that used in Fokker-Planck descriptions; as well as the algebraic relations among the 
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FIG. 6: Marginal distributions for the (a) first and (b) second gene in a two-gene network, summed from the joint distribution in 
Fig 5(b). Agreement is demonstrated among distributions obtained by the spectral method (solid), by direct iterative solution 
of the master equation (Eqn. 163; plus signs), and by standard Gillespie simulation [48] (circles). 

different cigenfunctions of the birth-death process itself, yielding an algebra of "ladder operators" similar to that used 
in quantum mechanics. 

While the literature in simulation of intrinsic noise has grown immensely since the introduction of time-varying 
Monte Carlo methods over thirty years ago [48], it is our hope that for model systems these calculations will help 
point out how the linearity of master equations may be exploited to yield either solutions, numerical methods, or 
analytical approximations which give more direct insight into the deign principles of stochastic regulatory networks. 



Appendix A: Deriving the master equation 



Here we derive the master equation in one dimension from the laws of probability. The probability of having n 
proteins at a time t + r (where r is small) is equal to the probability of having had n' proteins at time t multiplied 
by the (possibly time-dependent) probability of transitioning from n' to n in time r, summed over all n': 



V 



n(t + T) = J~]p n '(t)p n \ n i(t,T). 



(Al) 



The transition probability p n \ n i (t, r) has two contributions: (i) the probability of transitioning from state n' =/= n 
to state n, equal to a transition rate w nn '(t) times the transition time r, and (ii) the probability of starting and 
remaining in state n' = n, which we call 7r n (£): 



Pn\n'(t, T ) = (1 - <>rm')w nn '(t)T + 5 nn 'TT n (t) 

Applying the normalization condition ^2 n p n \n'(t> T ) = 1 to Eqn. A2, we obtain 

n n (t) = 1-r ^2 w n'n(t), 



(A2) 



(A3) 
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making Eqn. Al 



Pn(t + T) = ^Pn'(<) 



(1 - 5 nn >)w nn ,(t)T + 8 nn , 1 -T ^ U>n"n(*) 



n" ^n 



Ty^Pn>(t)w nn >(t) - Tp n (t)w nn + p n (t) - TP n (t) y~] W n n n {t) 

n' n"^n 



Pn{t) +T 



^2 Pn' (*)«W (*) - Pn (*) X! W ™'" (*) 



Taking the limit t — > 0, we obtain 



as in Eqn. 4. 






= X! P ™' Wn ™' (0 -Pn^W w '„ft), 



(A4) 
(A5) 

(A6) 
(A7) 



Appendix B: Summation of the birth— death master equation 

Here we sum the one-dimensional birth-death master equation over n against (i) n, to derive the mean dynamics, 
(ii) x n , to derive the equation of motion in generating function space, and (iii) \n), to find the equation of motion in 
operator notation. 

First we sum the master equation (Eqn. 9) over n against n to derive the dynamics of the mean (n) — ^ n np n : 



d t ^ n P« = ~9 X! n P« ~ X! n2 P™ + 9 X! n Vn-\ + ^2 n(n + l)p n +i- 



(Bl) 



The left-hand side (LHS) is the time derivative of the mean. On the right-hand side (RHS), recalling that n runs 
from to oo, we may define for the third term n' = n — 1 and for the fourth term n" = n + 1, giving 



d t (n) = -g J2 n Pn-J2 n2pn + 9 J2 ( n ' + 1 ) p "' + J2 ( n " _ W* 



(B2) 



fl' = -l 



The lower limits on all sums can be set to once more without changing the expression (in the third term we impose 
P—i = since protein number cannot be negative); simplifying gives 

d t (n) = -g ^ n Pn - X! n2pn + 9 X! n ' Pn ' + g X! Pn ' + X! n " 2 Pn" ~ ^ n " p ™" ( B3 ) 

n n n' n' n" n" 

= S^P„--^«V (B4) 

n' n" 

= 9-{n), (B5) 

as in Eqn. 8. 

Next we sum the master equation (Eqn. 9) over n against x n to derive the equation of motion of the generating 
function G(x,t) = J2 n p n {t)x n : 



d t ^p n x n = -g^2p n x n - ^np n x n + g^p„_ix" + ^(n+ l)p„+is 



(B6) 



The LHS and the first term on the RHS reduce directly to functions of G. The third and fourth terms on the RHS 
benefit from the same index shifts applied to obtain Eqn. B2, giving 



d t G= -gG - Y^ Pnnx n + g ^ p n >x n ' +1 + ^ n"p n 



„«"-! 



(B7) 



27 

We may now eliminate bare appearances of n in the second and fourth terms on the RHS using nx n = xd x x n and 
ra n_1 = d x x n respectively, giving 

d t G = -gG-Y / Pnxd x x n + gY,Pn'X n ' +1 + J2Pn"9 x x n " (B8) 

n n' n" 

= -gG - xd x G + gxG + d x G (B9) 

= -( x -l)(d x -g)G (BIO) 

as in Eqn. 23. 

Finally we sum the master equation (Eqn. 9) over n against \n) to derive the equation of motion of the generating 
function in operator notation \G) — ^2 n p n \n): 

dt^2p n \n) = -g^2p n \n) -^np n \n) + g^p n -\\n) + £(n + l)p„+i\n). (Bll) 

n n n n n 

Again we may reduce the LHS and the first term on the RHS directly, and we may apply index shifts to the third 
and fourth terms on the RHS: 

d t \G) = -g\G) -J2 n Pn\n)+gJ2Pn'\ n ' + 1) + $>' W - 1>- ( B12 ) 

n n' n" 

Now defining operators a + and a~ that raise and lower the protein number by 1 respectively, i.e. 

a+\n) = |n + l), (B13) 

cT\n) = n\n-l), (B14) 

Eqn. Bll can be written 

d t \G) = -g\G)-J2Pna + a-\n)+gJ2Pn'a + \n')+J2Pn"a-\n") (B15) 

n n ! n" 

= -g\G)-a+a-\G)+ga+\G)+a-\G) (B16) 

= -(a+ - l)(a~ - g)\G) (B17) 

as in Eqn. 50. 

Appendix C: Statistics of the birth— death distribution 

In this appendix we calculate the means of a birth-death distribution. In the simplest case, starting directly from 
the solution in Eqn. 16 we obtain: 



(„) = J2 nPn = £ n y -e-° = e -°gd g £ ^ = 2- ( C1 ) 

n n n 

Similarly we can calculate 

(n(n-l)) = Y, n ( n - l ^ 9 ^^ 9 9 2 d 2 g T,^^9 2 , (C2) 

n n 

which gives us the variance 

(n 2 )-(n) 2 = g 2 +g-g 2 =g. (C3) 

These results can also be obtained in the operator formalism, and the fact that for steady state \G) = \j = 0). For 
example: 

(n) = (x\a + a~\j = Q)\ x =\ = {x\a + oT \G)\ x= i = {x\a + aT £p„|n)| x= i = ^np n x n \ x=1 = £np„. (C4) 
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To use the number operator on the \j = 0) state we need to rewrite it in terms of a + = b + + 1 and aT = b~ + g. 
Acting to the left we obtain: 

(n) = (x\ (b+b- +g + b-+ gb+^j |0) x=1 = (x\g\0) x=1 + (x\g\l) x= i = g. (C5) 

We can also calculate the mean of the bursty process. Multiplying Eqn. 147 by n and summing over n, we use, 

oo oo 

n=0 n=0 

to obtain: 

g((n) + N)- g(n) + (n 2 ) - (n) - (n 2 ) = 0, (C7) 

which is solved by: 

(n) = gN. (C8) 

The second moment (n 2 ) = Yl n2 Pn is calculated the same way: 

(n 2 ) = \ (2g(n)N + (n) + gN 2 ) = (gN) 2 + ^ ± gN , (C9) 

which yields the variance as 

5n 2 = (n 2 )-(n) 2 = ^-(l + N). (CIO) 



Appendix D: Asymptotic distributions for large protein number 

Here we show that in the limit of large protein number n, the steady state of the birth-death process (the Poisson 
distribution, Eqn. 16) approaches a Gaussian distribution with mean and variance equal to g, the ratio of production 
to degradation rates. We then show that in the same limit, the steady state solution to the birth-death Fokker-Planck 
equation (Eqn. 84) also approaches a Gaussian distribution with mean and variance equal to g. 

The large-n limit of the Poisson distribution is most conveniently evaluated by first taking the log, which allows 
one to make use of Stirling's approximation, logn! w nlogn — n + log \2irn for large n: 

logp(n) rj — g + nlogg — nlogn + n — log v 2im. (Dl) 

The derivative 

d n logp{n) = log^-+0{l/n) (D2) 

n 

vanishes at the maximum 

n = g (D3) 

about which we Taylor expand to second order: 

logp(n) w logp(g) + -(n-g) 2 dl [log p(n)] g (D4) 

z (n — q) 2 , 

= -log y/^Tg- K y> . D5 

The result is a Gaussian distribution with mean and variance equal to g: 

p(n) = _^ e -(«-9) 2 /(2ff). (D6) 

y/2irg 
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The large-n limit of the Fokker-Planck distribution is similarly evaluated: the log and it's derivatives are 



N ( n\ 

logp(n) - log - + (40 - 1) log ( 1 + - , (D7) 

d n \ogp(n) = fc^-2, (D8) 

n + <? 



d*logp(n) = — * (D9) 



The first derivative vanishes at n = g — 1/2 ss g for large g, at which the second derivative evaluates to d\ [logp(n)] ( . 
— 1/g + l/(4g 2 ) w —1/g. Taylor expanding to second order and exponentiating then give 



p(n) = — ^ e -("-9) 2 /(2ff) 7 (dio) 

V27Tff 



where TV is eliminated by normalization. 



Appendix E: Orthonormality and the inner product 

Here we show that interpreting the generating function as a Fourier transform motivates a particular choice of inner 
product and conjugate state in the protein number basis. 

We have defined the generating function in terms of a continuous variable x as G{x) = ^2 n Pn% n - We could 
equivalently write x = e lK and define 

G(K) = ^PnJ™, (El) 

n 

which makes clear that the generating function is simply the Fourier transform of the probability distribution in protein 
number n. In the state notation commonly used in quantum mechanics [42], Eqn. El is written (k\G) = ^2 n Pn{n\n) , 
where (n\G) = G(k) and 

(«|n) = e inK . (E2) 

Eqn. E2 is the representation of \n) in n space. In order to compute projections of \n) onto other states or itself 
using this representation, we must define a conjugate state (ti\k) and a consistent inner product. With complex 
exponentials, it is common to make the conjugate state the complex conjugate, 

(n|«) = («|n)* = e~ inK , (E3) 

and the inner product the integral 

(f\h)=j -(.m( K \h) 7 (E4) 

for any / and h. Given the definition of conjugate state, this choice of inner product ensures the orthonormality of 
the \n) states: 

p-(n\K)(n\ri) = / ^ n "")" = S nn ,. (E5) 

o 27r Jo 27T 

We may now reinterpret these definitions in terms of our original variable x. Since x = e lK we have dx — ie tK dK = 
ixdn, and in the orthonormality condition (Eqn. E5) the integration from to 2tt along the real k line becomes contour 
integration along the unit circle in the complex x plane: 

(n\n') = r *e-» fc e'"' fc - / ^ x - n x n ' = / p,-^x n ' = 6 nn ,. (E6) 



/ 2tt J 2-kix J 2iri x n+1 

Since we have defined 

(x\n) = x n , (E7) 
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Eqn. E6 suggests the definition of conjugate state 

("I*) = Jtt ( E8 ) 

and inner product 

(f\h) = f ^(f\x)(x\h) (E9) 

in x space. 

With these definitions, we will find Cauchy's theorem, 

^(^1 = ^^-^ (E10) 

(with the convention 9(0) = 1 for the Heaviside function) very useful in evaluating projections. For example, we may 
immediately use it to confirm the orthonormality condition, 

/dr 1 / 1 r n 

2^^n x " = d a " r \ x= o e{n) = <5n -'' (E11) 

and in Appendix H we use it to compute analytic expressions for the projections between protein number states \n) 
and birth-death eigenstatcs \j). 

Appendix F: Properties of the raising and lowering operators 

Just as in the operator treatment of the quantum harmonic oscillator [42] , in this paper we have defined operators 
a + and a~ that act on \n) states by raising and lowering the protein number n by 1 respectively, i.e. 

a+\n) = |n+l), (Fl) 

ar\n) = n\n- 1) (F2) 

(note, however, that the prefactors here, 1 and n, are different than those conventionally used for the harmonic 
oscillator, \pa + 1 and ^/n, respectively). In this appendix we derive the actions of these operators when acting to 
the left, as well as their commutation relation. 

The actions of a + and a~ to the left can be found by projecting onto Eqns. Fl and F2 the conjugate state (n'\: 

(n'\d+\n) = (n'\n + 1) = 6 n >, n+1 = S n >-i, n - (n' - l|n), (F3) 

(n'\a-\n) = (n'\n\n - 1) = rwJ„/, n _i = (n' + l)6 n >+i, n = (n' + l)<n' + l|n). (F4) 

The first and last expressions in each case imply 

(n'\a + = (n'-l|, (F5) 

(n'\a- = (n' + l)(n' + l|, (F6) 

as in Eqns. 56-57. 

The commutation relation between a~ and a + is defined as [aT , d + ] = aT a + — a + aT . It is evaluated by considering 
its action on a state \n): 

[a~, a + ] \n) = a~ a + \n) — a + a~\n) = a~\n + 1) — a + n\n — 1) = (n + l)|n) — n\n) = \n). (F7) 

The first and last expressions imply 

[a-,o+]=l, (F8) 

just as with the quantum harmonic oscillator. 
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Appendix G: Eigenvalues and eigenfunctions in the operator representation 

Here we derive the eigenfunctions \Xj) (also called eigenstates in the operator representation) and eigenvalues Xj of 



the operator C = b + b , for which 



C\X j ) = X j \X j ), 



(Gl) 



as well as the actions of the individual operators b + and b on the eigenstates. We will find that the commutation 
relation 



b-y 



1 



and the existence of the steady state solution (x\G) — e qlyX *), for which 

t\G) = 0, 



(G2) 



(G3) 



are all that are necessary to completely define the eigenvalues and eigenstates of C The treatment here finds many 
parallels with the derivation of the eigenvalue spectrum of the quantum harmonic oscillator [42] . 
First it is useful to compute the commutation relation of C with each of its components b + and b~ : 



~b+,t 

b~,l 


— 


b + ,b + b~ 

h-,b+b- 


= b+ 


b + ,b- 
h-,i>- 


+ 
+ 


'b+,b+ 

b~,b + 



b- = -F 



(G4) 
(G5) 



(here we have used Eqn. G2 and the properties [/, gh] — g[f, h] + [f,g]h and [/, /] = for any /, g, and h). We now 
consider a particular eigenvalue A, and evaluate the action of Lb + on |A): 



b + X 



Cb+\X) = b+C\X) - 
The equality of the first and last expressions, 

£[b + \x)) = (x + i)(is\x) 



\X)=b+X\X)+b+\X) = (X + l)'b+\X). 



(G6) 



(G7) 



reveals two results: (i) the existence of a state with eigenvalue A implies the existence of a state with eigenvalue 
A + 1, and (ii) the state with eigenvalue A + 1 is proportional to b + \X). By induction the first result means that the 
eigenvalues are spaced by 1; moreover, the existence of a state with eigenvalue zero (the steady state solution; Eqn. 
G3) anchors the eigenvalues to the integers, i.e. 



x 3 = J 
for integer j. The second result can now be written explicitly 

h + \j) = \j + i), 



(G8) 



(G9) 



where we are free by normalization to set the proportionality constant to 1. Eqn. G9 demonstrates that b + is a raising 
operator for the eigenstates \j). 

Evaluating the action of OoT on cigenstate \j) similarly reveals that b~\j) is proportional to \j — 1) (cf. Eqn. G6); 
consistency with the eigenvalue equation (Eqn. Gl) requires that the proportionality constant be j, giving 



b-\j)=j\j-l). 



(G10) 



Eqn. G10 demonstrates that b is a lowering operator for the eigenstates \j). The operators b + and b raise and 
lower \j) as a + and aT do |n); therefore they act to the left as (cf. Eqns. F3-F4) 



U\b + = (i-l|, 

U\b- = (j + l)0' + l|. 



(Gil) 
(G12) 
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Eqn. G10 imposes a floor on the eigenvalue spectrum, since fo |0) = 0, and no lower eigenstates can be generated. 
Therefore the eigenvalues are limited to nonnegative integers: 



j €{0,1,2,3,...}. 
Eqn. G9 describes how any state can be obtained from the j = state: 



(G13) 



(G14) 



Recalling that b + = a + — 1 and that |0) is the steady state solution, Eqn. G14 can be used to derive the representation 
of the eigenfunctions in x space, (x\j). Projecting (x\ onto Eqn. G14 and recalling that a + corresponds to x (Eqn. 
54) and that (x|0) = e s ( x_1) (Eqn. 24), we obtain 



(x\j) = (s-l) J 'e»< a, - 1 >. 



(G15) 



We may now define the conjugate state (j\x) such that the \j) are orthonormal under the inner product defined in 
Eqn. E9: 

Sir = 01/> = / |^0W W> - / ^i<3\xKx iy'e^-v = I ^y j 'fi(y), (Gie) 



where y = x—1 and fj(x—l) = e 9 ^ x l > (j\x). Eqn. G16 is equivalent to Eqn. E6, from which we identify fj(y) = 1/y- 



j'+i 



and thus 



m = 



e -g(x-l) 

(x - l)J'+i ■ 



(G17) 



We have shown that the operator C = b + b has nonnegative integer eigenvalues j and eigenfunctions (in x space) 
given by Eqns. G15 and G17. 



Appendix H: Overlaps between protein number states and eigenstates 

Here we describe two methods for computing the overlaps (n\j) and (j\ri) between the protein number states \n) 
and the eigenstates \j): by contour integration and by recursive updating. 

The first method evaluates the overlaps using the inner product defined in Eqn. E9 and Cauchy's theorem (Eqn. 
E10). Recalling the representations in x space of the protein number states and eigenstates (Eqns. 46, 48, and 65-66), 
the first overlap becomes 



(n\j) 



dx 
2vri 



'n\x){x\j) 



dx e 9 ^ x -^(x-l) j 



2m 



„n+l 



— d n 
n! x 



p9(x-1) 



(x - If 



x=Q 



Repeated derivatives of a product follow a binomial expansion, i.e. 

pS(^-I) 



1 n 1 

Hi) = ^E^fT)! 9 "" 



E 



i 



e\(n-£)\ 



[9 n - e e-°] 



where we define 



1=0 

= (-l) j e- 9 g n Mni, 



min(n,j) 

S"j — 2^ 

£=Q 



(j-e) 



M^-^Uo 



i-iy- £ 9(j-£) 



£\(n-£)\(j-£)\(-gy 



Following similar steps, the conjugate overlap evaluates to 

(j\n) = n\(-g)°£, n j, 



(HI) 

(H2) 

(H3) 

(H4) 

(H5) 
(H6) 
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with £ n j as in Eqn. H5. For the special case j = 0, Eqns. H4 and H6 reduce to Eqns. 77-78. 

It is more computationally efficient to compute the overlaps recursively using rules that can be derived from the 
raising and lowering operations (Eqns. 52-53, 56-57, and 67-70). For example, using the raising operators, 

(n\j + 1) = (n\b+\j) = (n\(a+ - l)|j) = (n - l\j) - (n\j), (H7) 

which can be initialized using (n|0) = e~ 9 g n jn\ (Eqn. H4) and updated recursively in j. Eqn. H7 makes clear that in 
n space the (j + l)th mode is simply the (negative of the) discrete derivative of the jth mode. The lowering operators 
give an alternative update rule, 

(n + l)(n + l\j) = (n\a-\j) = <n| (S~ + g)\j) = j(n\j - 1) + g(n\j), (H8) 

which can be initialized using (0\j) = (— l) J e~ 9 (Eqn. H4) and updated recursively in n. One may similarly derive 
recursion relations for (j\n), i.e. 

(j\n + l) = (i-l|n) + 01n), (H9) 

(j + l)(j + l\n) = n (j\n-l)-g(j\n), (H10) 

initialized with (j|0) = {—g) J /j\ or (0|n) = 1 respectively (Eqn. H6) and updated recursively in n or j respectively. 
Two-term recursion relations can be similarly derived from the full operator b + b~ [29]. 

Appendix I: The equivalence of the Fokker— Planck and Langevin descriptions 

In this appendix we start from the Langevin equation defined in Eqns. 93 and 94 and derive the Fokker-Planck 
equation in Eqn. 81. As described in Sec. II D the observed trajectory is one realization r, q (i) of the random stochastic 
process with Gaussian noise as presented in Eqn. 92: 

p n (t)= jVr 1 5{n-r r} {t))P[ri\:={5{n-r r ,{t))). (II) 

Consider the evolution of the probability distribution: 

p n {t + At) - Pn (t) = (6(n- r v {t + At)) -5{n- r v {t))), (12) 

and expand the increment in the trajectory as r(t + At) = r(t) + Ar(t). We can now Taylor expand the difference in 
delta functions: 

Pn (t + At)- Pn (t) = ((-Ar(t))6'(n-r 11 (t)) + ^(-Ar(t)) 2 S"(n-r v (t))), (13) 

= -d n (Ar(t)5 (n r„(*))) + l -d 2 n ({Ar(t)f5 (n r„(t))>, (14) 

where the primes denote derivatives in n. Using the Langevin equation to calculate the increments in the trajectories: 

Ar(t) = v [r(t)] At + rj(t)At (15) 

we obtain: 

(Ar(t)S (n - r n (t))) = (v [r(t)\ Atd (n - r„(t))> + (t)(t)At5 (n - r„(t))) (16) 

= v [r(t)] At(5 (n - r„(t))) = v (n])At Pn (t) (17) 

((Ar(t)f5 (n r v (t))) = (At) 2 ((v [r(t)]fS (n r,(t))) + 2(v [r(t)} V (t)S (n r„(t))) + {r,\t)5 (n r,(t))) (18) 

= (At) 2 (v(n)) 2 Pn (t) + (r 1 2 (t)) Pn (t), (19) 

where we have used (assuming a discretized process): 

( V (t)5 (n - r,(t))) = (v(t))(S (n - r„(t))) = 0. (110) 

We identify 

D{n) = ^{r, 2 {t)). (Ill) 
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Putting together all elements and keeping only leading terms in At: 

Pn (t + AQ- Pn (t) = ^ [v {n])Atpn{t)] + Ql [D{n)pn{t)] . (I12) 

Taking the limit of At — > we recover the Fokker-Planck equation: 

dtPn{t) = -dn [v (n])At Pn (t)} + d 2 n [D{n)p n {t)\ . (113) 

Appendix J: Derivation of the Hill function 

The Hill function can be derived as an effective production rate for an autoregulating gene with two production 
states (Sec. IV A 4). The gene is found either in an inactive state (— ), in which the production rate is a constant g_, 
or in an active state (+), in which the production rate oj+n h depends on the number of proteins n, which incorporates 
the autoregulation; h describes the cooperativity, with h > corresponding to activation and h < corresponding to 
repression. The probability p^ of the gene being in a given state + or — and there being n proteins evolves in time 
according to the master equation in Eqn. 117 with £l zz i as in Eqn. 139; in steady state: 

= -£ ± p±±uj + n h p-TU-p+ 7 (Jl) 

where C^ describe simple birth-death terms with constant production rates g± in each of the two states. We define 
moments of the master equation as 

^ = $>+, (J2) 

n 

l^fif = J2 P n ni fOT ( ^ 1 - ( J3 ) 

n 

with 

7T+ + 7T- = 1 (J4) 

by normalization. 

Summing Eqn. Jl (top signs) over n (and recalling that the birth-death terms sum to zero) gives 

= w + Y,n h P~ -w_5>+ (J5) 

n n 

= UJ + 1T~ fl^ — 0J-7T + (J6) 

which, with Eqn. J4, becomes an expression for 7r + , the probability of being in the active state: 

7T+ = - f h . (J7) 

We solve for fxT using two approximations. The first is that higher moments can be be decoupled, i.e. 

Vh+i ~ HVi for h>l, (J8) 

which implies that 

% « (Mr)' 1 - (J9) 

This approximation allows one to simplify the mean equation in the + state (obtained by summing the Eqn. Jl, top 
signs, against n over n): 

= 7T+5+ - J2 n Pt + W + J2 ^^Pn - U -J2 n Pn ( J10 ) 

n n n 

= TT + g+ — 7T + /i+ + CJ + 7r _ /l^ +1 - CJ_7T + ^ (JH) 

« n + g + — n + fi~l + cj + 7t~ l\^i — uJ-ir + ^ (J12) 

= 7T + .9+ - 7r + A^ + w_7T+ (/z7/ - f4) , (J13) 
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where the last step uses Eqn. J6. The second approximation is that transitions between states are fast, i.e. uj + ~ 
cj_ 3> 1. This approximation allows us to neglect the first two terms of Eqn. J13 compared to the third term, which 
implies that /i^ « /i^ . Summing Eqn. J3 over ± for £ = 1 then gives 

tt + i±1 + -K-lXy « (n + + 7r~) [i{ = fi^ = y^p„n = n, (J14) 

n 

which shows that fj,^ approximates the mean of the distribution. Therefore with Eqn. J9 the probability that the 
gene is active (Eqn. J7) can be written 

(J15) 



n h +K' 



with equilibrium constant K = lu + /lu-. The effective production rate is the sum of the production rates in each state 
times the corresponding probabilities of being in each state: 

g(n) = g-Tr~+g + ir + (J16) 

(J17) 



which is the Hill function (Eqn. 108). 



, (l ~ nh ^ 


n h 


g-K + g+n h 


1 ' 9+ n* + K 



, ,- ■ (J18) 

n h + K ' 



Appendix K: Limiting case of the two-state gene 

Here we show that the spectral solution to transcriptional bursting (Eqn. 137) reduces to the hypergcomctric form 
(Eqn. 127) in the limit of Z — 2 production states. We also derive the slightly simpler expression for the special case 
of zero production in the inactive state. 

In the case of two states (Eqns. 120-121; z = ±), Eqn. 137 becomes 

jGf + u T Gf - ^Gj =u ± Y: GJ, ^9f, (Kl) 

where A = A^ = —A |_. Initializing with Gq = uj±/(lj + + w_) and computing the first few terms reveals the 

pattern 

G i = - — rr. ^— „,--i ,... . : rr r ( K2 ) 



u± (TA) J T(j + ufcp) T{lu+ + w_ + 1) 



W_|_ + w_ 



j\ T(cu T ) T(j + LU+ + LO_ + 1) ' 



(K3) 



where in the second line the products are written in terms of the Gamma function. Writing the total generating 
function \G) = Yl± \G±) in position space recovers the hypergeometric form (Eqn. 127): 

G(x) = "£{x\G ± ) (K4) 

± 

= EEw±>o±ig ± > (K5) 

± 3 

= ^^{x-iye g ^ x -^Gf (K6) 

± 3 

= y _^ e^ x - 1 1^[Lu T ,Lo + +Lo_ + l; T A(x-l)}, (K7) 
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where 

$\a8-u]=y r{j + a) m U - (K8) 

9[a,p,u\ ^ r(a) r(?+/3)? , {**) 

is the confluent hypergcomctric function of the first kind. 
In the limit g_ = 0, Eqn. K7 reads 

G(x) = e"$[w_,u; + +w_ + 1; -u] H = ^\ui + ,ui + +u>- + l;u], (K9) 

where u = ,g+(.x — 1). Using the fact that [49] 

e u $[a,/3;-w] = $[/?- a, /3;u], (K10) 

Eqn. K9 can be written 

G(x) = — — $[w+ + l,w++w_ + l;d H — $[w + ,cj + +w_ + l;w], (Kll) 

or, noting Eqn. K8 and the fact that T(s + 1) = sT(s) for any s, 

G(x) = W ^+ T U+ UJ + + 1 ) u- r(j + o; + ) \ I> + +^_ + l) v? 

{> ^\uo+ + uj- r(w+ + i) uj + + uj- r(w+) y r(i + w+ + cj_ + i) j!' l j 

V ( UJ+ ^' + ^+) r ( J+u; +) | UJ - T U + UJ +) \ (uj++uj-)T(uj + +uj-) u> 

/ -^\u; + +u;_ lo + T(lu + ) oj + + cj_ r(w + ) / (j + lu + + LU-)T(j + uj + + u>_) j! ' 

\p r(j + g+) r(^ + + ^>_) M J 

^ r( W+ ) r(j + W+ +a;_)i! l j 

= $[w+,w++w_;w], (K15) 

as in Eqn. 128. The marginal p„ is obtained by p„ = d™[G(x)]o/n\; using Eqn. K15 and the derivative of the confluent 
hypergcomctric function, 

d:${a,{3;u} = r( " + a) r( f $[a + n,p + n;u], (K16) 

r(a) r(n + p) 



one obtains 



as in Eqn. 129. 



g r LT(n + uj + ) T(uj + +lu) , , , 

*» = l IV S r( m ,' M u++n,u++u-+n;-9+], K17 
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